In [48]:
# Create a list of AnnData objects (one for each sample)
adatas = []
for sample, sample_meta in samples.items():
print("PROCESSING...")
# Load the associated transcriptomics data
gex_file = glob(f"{data_dir}*{sample}*.h5")[0]
adata = sc.read_10x_h5(gex_file)
# Load the TCR data
tcr_file = glob(f"{data_dir}*{sample}*airr*")[0]
adata_tcr = ir.io.read_airr(tcr_file)
print(sample, " cells with tcr:", adata_tcr.shape)
print(sample, " cells with GEX:", adata.shape)
#Preprocessing
adata.var_names_make_unique()
sc.pl.highest_expr_genes(adata, n_top=20, save="_"+sample+"_highestGenes."+imageFormat)
#basic filtering
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
adata.var['mt'] = adata.var_names.str.startswith('MT-') # annotate the group of mitochondrial genes as 'mt'
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],
jitter=0.4, multi_panel=True, save="_"+sample+"_QCmetrics."+imageFormat)
#Remove cells that have too many mitochondrial genes expressed or too many total counts:
sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt', save="_"+sample+"_UMI_mt."+imageFormat)
sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts', save="_"+sample+"_UMI_feat."+imageFormat)
#Actually do the filtering by slicing the AnnData object
adata = adata[adata.obs.n_genes_by_counts > 300, :]
adata = adata[adata.obs.total_counts > 500, :]
adata = adata[adata.obs.total_counts < 10000, :]
adata = adata[adata.obs.pct_counts_mt < 10, :]
print("After filtering", sample, " cells with GEX:", adata.shape)
ir.pp.merge_with_ir(adata, adata_tcr)
adata.obs["sample"] = sample
adata.obs["group"] = sample_meta["group"]
# concatenation only works with unique gene names
adata.var_names_make_unique()
adatas.append(adata)
PROCESSING... BP_01 cells with tcr: (837, 0) BP_01 cells with GEX: (1348, 36601)
After filtering BP_01 cells with GEX: (1293, 13779) PROCESSING... BP_02 cells with tcr: (530, 0) BP_02 cells with GEX: (930, 36601)
After filtering BP_02 cells with GEX: (916, 12867) PROCESSING... BP_03 cells with tcr: (222, 0) BP_03 cells with GEX: (410, 36601)
After filtering BP_03 cells with GEX: (400, 11023)
